Ce document montre lāessentiel des manipulations Ć effectuer avec R pour rĆ©pondre aux questions du TP1 dāAnalyse UnivariĆ©e. Il suppose une connaissance prĆ©alable de la syntaxe R, de prĆ©fĆ©rence Ć lāaide de lāIDE RStudio et les bases du package sf (tutoriel ici) qui implĆ©mente la norme SF comme dans PostGIS.
La connaissance du package dplyr et de sa notation à base de pipe (%>%) sera utile pour écrire des chaînes de traitements plus concises que celles listées ici.
Nous avons rĆ©cupĆ©rĆ© le fichier des donnĆ©es spatialisĆ© sous le format GeoJSON. Dāautres formats sont disponibles, notamment via lāAPI du site. Nous allons utiliser le package sf pour manipuler les donnĆ©es spatiales. Pour lire le GeoJSON , nous avons besoin dāinstaller un package spĆ©cifique : geojsonsf, qui se charge de la conversion entre des objets geographiques au format GeoJSON et leur reprĆ©sentation au format sf.
Lāinstallation et le chargement des packages au moyen des fonctions install.packages("nom_du_package") et library(nom_du_package)
install.packages("geojsonsf")
library(geojsonsf)
install.packages("sf")
library(sf)
Nous installons et chargeons aussi les packages dplyr (manipulation aisƩe de donnƩes) et cartography (cartographie correcte des objets spatiaux).
install.packages("dplyr")
library(dplyr)
install.packages("cartography")
library(cartography)
Le fichier geoJSON est chargĆ© Ć lāaide de la fonction geojson_sf().
N.B. la fonction geojson_sf requiert le chemin dāaccĆØs au fichier, fourni sous sa forme relative ou absolue. Le repertoire de travail de R est dĆ©fini par la commande setwd("/chemin/vers/le/fichier"). Pour connaĆ®tre le repertoire de travail courant, utiliser la fonction getwd()
N.B. le nombre de lignes et les valeurs de données qui apparaissent dans ce support peuvent varier, les données étant mises à jour régulièrement sur le site opendata.paris.fr
arbres <- geojson_sf("les-arbres.geojson")
names(arbres) # affiche le nom des colonnes (variables)
## [1] "remarquable" "circonferenceencm" "stadedeveloppement"
## [4] "genre" "idbase" "arrondissement"
## [7] "idemplacement" "geo_point_2d" "geometry"
## [10] "adresse" "libellefrancais" "complementadresse"
## [13] "domanialite" "typeemplacement" "hauteurenm"
## [16] "varieteoucultivar" "espece"
Le chargement du fichier peut ĆŖtre un peu long: il contient 203818 observations le jour de la rĆ©daction de ce support. Pour connaitre la structure de lāobjet , nous utilisons la fonction str() qui nous donne un aperƧu du type des colonnes du dataframe que nous manipulons. Cāest un dataframe particulier puisquāil est Ć©galement de la classe sf : lāun de ses colonnes est dĆ©diĆ© Ć la description de sa gĆ©omĆ©trie (ici, des points 2D)
La fonction head() permet dāobserver les \(n\) premiĆØres lignes dāun dataframe. Ici la fonction paged_table est utilisĆ©e pour un affichage dynamique plus pratique du tableau dans ce support HTML.
paged_table(head(arbres, 10))
La librairie sf surcharge la fonction dāaffichage par dĆ©faut de R , plot, pour afficher la gĆ©omĆ©trie des donnĆ©es comme un objet gĆ©ographique et non comme un nuage de points ou une courbe. Auparavant , nous devons fixer le systĆØme de coordonnĆ©es de rĆ©fĆ©rences de lāobjet, pour que les donnĆ©es soient correctement projetĆ©es Ć lāaffichage.
La fonction st_crs() appliquĆ©e sur un objet spatial retourne son CRS sāil est dĆ©fini , ou permet de fixer avec lāopĆ©rateur dāaffectation <-. Nous allons utiliser Lambert 93 (EPSG 2154). LāopĆ©ration se fait en deux temps : transformation selon le CRS puis affectation de lāattribut.
st_crs(arbres)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
arbres <- st_transform(arbres,2154)
arbres <- st_set_crs(arbres, 2154)
st_crs(arbres) # nouveau CRS de l'objet aprĆØs transormation
## Coordinate Reference System:
## EPSG: 2154
## proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
Nous pouvons maintenant afficher les donnƩes avec la fonction plot
plot(arbres)
## Warning: plotting the first 9 out of 16 attributes; use max.plot = 16 to
## plot all
## Warning in classInt::classIntervals(na.omit(values), min(nbreaks, n.unq), :
## N is large, and some styles will run very slowly; sampling imposed
Cāest trĆØs long ! CelĆ est dĆ» au fait que la fonction plot affiche autant de cartes que de variables (colonnes) de lāobjet spatial, jusquāĆ un maximum de 9 colonnes apr dĆ©faut.
Pour le moment nous nāallons que tracer la gĆ©omĆ©trie des donnĆ©es : lāattribut geometry de lāobjet
plot(arbres$geometry, cex = 0.01)
Tout autre attribut peut être choisi , la fonction de dessin de R se charge de trouver une échelle de couleur en fonction des valeurs que prend la varaible.
plot(arbres["arrondissement"], cex = 0.01, graticule=T)
La sequence de commande suivantes transforme lāattributs arrondissement en facteurs et filtre les donnĆ©es de faƧon Ć ne conserver que les arrondissements de Paris intra muros.
Nous utilisons les fonctions:
as.factor() pour transformer une colonne en facteurs , i.e.Ā des variables qui ne peuvent prendre quāun certain nombre de modalitĆ©slevels() qui donnent les modalitĆ©s que peut prendre une variable factoriellefilter() ,fonction du pakcage dplyr qui ne conservent que les Ć©lĆ©ments dāun tableau correspondant Ć un certain prĆ©dicat boolĆ©en passĆ© en argument%in% (appartenance) et lāopĆ©rateur boolĆ©en ! (NOT)as.data.frame() qui convertit un objet en dataframe (lorsque cāest possible)droplevels , qui supprime les modalitĆ©s inutilisĆ©es dāun facteurarbres$arrondissement <- as.factor(arbres$arrondissement)
levels(arbres$arrondissement) # modalitƩs de la variable
## [1] "BOIS DE BOULOGNE" "BOIS DE VINCENNES" "HAUTS-DE-SEINE"
## [4] "PARIS 10E ARRDT" "PARIS 11E ARRDT" "PARIS 12E ARRDT"
## [7] "PARIS 13E ARRDT" "PARIS 14E ARRDT" "PARIS 15E ARRDT"
## [10] "PARIS 16E ARRDT" "PARIS 17E ARRDT" "PARIS 18E ARRDT"
## [13] "PARIS 19E ARRDT" "PARIS 1ER ARRDT" "PARIS 20E ARRDT"
## [16] "PARIS 2E ARRDT" "PARIS 3E ARRDT" "PARIS 4E ARRDT"
## [19] "PARIS 5E ARRDT" "PARIS 6E ARRDT" "PARIS 7E ARRDT"
## [22] "PARIS 8E ARRDT" "PARIS 9E ARRDT" "SEINE-SAINT-DENIS"
## [25] "VAL-DE-MARNE"
arbres_intramuros <- filter(arbres, !(arrondissement %in% c("BOIS DE BOULOGNE", "BOIS DE VINCENNES", "HAUTS-DE-SEINE", "SEINE-SAINT-DENIS", "VAL-DE-MARNE")))
# on retire les valeurs modales inusitƩes de la variable
arbres_intramuros$arrondissement <- as.character(arbres_intramuros$arrondissement)
arbres_intramuros$arrondissement <- as.factor(arbres_intramuros$arrondissement)
levels(arbres_intramuros$arrondissement)
## [1] "PARIS 10E ARRDT" "PARIS 11E ARRDT" "PARIS 12E ARRDT"
## [4] "PARIS 13E ARRDT" "PARIS 14E ARRDT" "PARIS 15E ARRDT"
## [7] "PARIS 16E ARRDT" "PARIS 17E ARRDT" "PARIS 18E ARRDT"
## [10] "PARIS 19E ARRDT" "PARIS 1ER ARRDT" "PARIS 20E ARRDT"
## [13] "PARIS 2E ARRDT" "PARIS 3E ARRDT" "PARIS 4E ARRDT"
## [16] "PARIS 5E ARRDT" "PARIS 6E ARRDT" "PARIS 7E ARRDT"
## [19] "PARIS 8E ARRDT" "PARIS 9E ARRDT"
On peut maintenant tracer les positions des arbres de Paris intra-muros, après avoir vérifié le CRS de cette couche par la fonction st_crs() (pour superposer plus tard les deux couches si celles-ci )
st_crs(arbres_intramuros)
## Coordinate Reference System:
## EPSG: 2154
## proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
plot(arbres_intramuros["arrondissement"], cex=0.1)
Nous allons maintenant charger les contours des quartiers de Paris, disponibles sur (https://opendata.paris.fr/explore/dataset/quartier_paris/information/). Chaque arrondissement contient 4 quartiers, on a donc 80 unités spatiales à considérer.
quartiers <-(read_sf("quartier_paris.shp"))
quartiers <- st_transform(quartiers, 2154)
quartiers <- st_set_crs(quartiers,2154)
plot(quartiers$geometry)
plot(arbres_intramuros["arrondissement"], add=TRUE, alpha=0.5, cex=0.1 )
Pour tracer un objet spatial , il suffit dāappliquer la fonction plot() sur cet objet. Cette fonction peut Ć©galement colorer la gĆ©omĆ©trie de lāobjet en fonction dāune variable.
Nous utilisons les fonctions count()du package dplyr, pour compter les nombre dāarbre par genre. (cāest lāĆ©quivalent dāun group_by, suivi dāune aggregation de comptage : count) Nous utilisons la fonction top_n du package dplyr pour selectionner les 6 genres les plus reprĆ©sentĆ©s. Nous utilisons la fonction filter et lāoperateur %in% pour ne conserver que les arbres de ces six genres.
count_by_genre <- count(arbres_intramuros,genre, libellefrancais, sort = T,name = "nb_indiv")
head(count_by_genre, 10)
## Simple feature collection with 10 features and 3 fields
## geometry type: MULTIPOINT
## dimension: XY
## bbox: xmin: 644413.8 ymin: 6857489 xmax: 657170.6 ymax: 6867050
## epsg (SRID): 2154
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 10 x 4
## genre libellefrancais nb_indiv geometry
## <chr> <chr> <int> <MULTIPOINT [m]>
## 1 Plata⦠Platane 37856 (645035 6861100, 645068.5 6861382, 645ā¦
## 2 Aescu⦠Marronnier 20051 (644413.8 6861116, 644420.5 6861114, 6ā¦
## 3 Tilia Tilleul 16778 (645059.1 6860934, 645070.1 6861108, 6ā¦
## 4 Acer Erable 12369 (645032.2 6860884, 645037.9 6860896, 6ā¦
## 5 Sopho⦠Sophora 10708 (645069.1 6861372, 645088.8 6860152, 6ā¦
## 6 Pinus Pin 3751 (645037.9 6860889, 645057.3 6861095, 6ā¦
## 7 Celtis Micocoulier 3628 (645102.9 6860998, 645167.5 6861102, 6ā¦
## 8 Fraxi⦠FrĆŖne 3580 (645068.1 6860577, 645074.4 6861408, 6ā¦
## 9 Prunus Cerisier Ć fleu⦠3359 (645037.9 6860884, 645056.4 6860576, 6ā¦
## 10 Pyrus Poirier Ć fleurs 3230 (645093.6 6861056, 645103.2 6861005, 6ā¦
six_genres <- top_n(count_by_genre,6,nb_indiv)
head(six_genres)
## Simple feature collection with 6 features and 3 fields
## geometry type: MULTIPOINT
## dimension: XY
## bbox: xmin: 644413.8 ymin: 6857489 xmax: 657126.1 ymax: 6867050
## epsg (SRID): 2154
## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 6 x 4
## genre libellefrancais nb_indiv geometry
## <chr> <chr> <int> <MULTIPOINT [m]>
## 1 Plata⦠Platane 37856 (645035 6861100, 645068.5 6861382, 64509ā¦
## 2 Aescu⦠Marronnier 20051 (644413.8 6861116, 644420.5 6861114, 644ā¦
## 3 Tilia Tilleul 16778 (645059.1 6860934, 645070.1 6861108, 645ā¦
## 4 Acer Erable 12369 (645032.2 6860884, 645037.9 6860896, 645ā¦
## 5 Sopho⦠Sophora 10708 (645069.1 6861372, 645088.8 6860152, 645ā¦
## 6 Pinus Pin 3751 (645037.9 6860889, 645057.3 6861095, 645ā¦
Lāimmense avantage du package sf : les gĆ©omĆ©tries ont Ć©tĆ© conservĆ©es lors des manipulation de regrouppement , filtrage etc. Il nāest donc pas nĆ©cessaire de refaire les opĆ©rations sur les gĆ©omĆ©tries , qui ont āsuiviā les rĆ©sultats lors des opĆ©rations count et top_n .
plot(six_genres["libellefrancais"],
key.pos=1, cex=0.2,
key.length = 0.9,
main="Les six genres d'arbres les plus représentés à Paris")
Que peut-on dire de lāimplantation de ces 6 genres ?
Afficher cette variable est immĆ©diat. La lĆ©gende est perfectible pour le moment; pour faire une vraie carte , voir Ć la fin du support lāutilisation du package cartography.
plot(arbres_intramuros["domanialite"],
key.pos=4, cex=0.2,
key.length = 1,
key.width = lcm(4.5),
main="DomanialitƩ des arbres de Paris")
On rĆ©alise une jointure (spatiale) entre quartiers et arbres_intramuros, pour inclure les informations de quartiers aux arbres Ć lāaide du prĆ©dicat spatial st_within
arbres_in_quartiers <- st_join(arbres_intramuros,quartiers, join=st_within)
nrow(arbres_in_quartiers)
## [1] 163367
On obtient alors un nouveau dataframe spatial, contenant les 163367 arbres situƩs dans les quartiers intramuros, et augmentƩs des variables issues du dataframe quartier.
On peut alors compter le nombre dāarbres par quartiers avec la fonction table qui compte (entre autres) le nombre dāindividus (ligne) par valeurs distinctes dāune de ses variables. On peut Ć©tendre le nombre de variables pour obtenir des tables de contingence.
nb_arbres_by_quartier <- table(arbres_in_quartiers$c_qu)
Pour ajouter cette information au dataframe quartiers, lāune des possibilitĆ© est de trier le dataframe par son code de quartier c_qu puis dāajouter la colonne nb_arbres_by_quartiers obtenue prĆ©cĆ©demment.
on peut par exemple utiliser la fonction arrange du package dplyr
quartiers <- arrange(quartiers,c_qu)
quartiers$nb_arbres <- nb_arbres_by_quartier
Avant de cartographier cette valeur, regardons sa distribution Ć lāaide dāun histogramme pour dĆ©terminer si un traitemment particulier doit ĆŖtre envisagĆ© (en cas de distribution particuliĆØrement biscornue)
hist(quartiers$nb_arbres,breaks = 10)
plot(quartiers["nb_arbres"], main=NULL, key.pos = 4)
title("Nombre d'arbres par quartier administrif de Paris",sub = "ReprƩsentation sƩmiologiquement discutable" )
Pour rĆ©aliser une sorte de raster, nous allons crĆ©er une grille qui couvre la zone dāĆ©tude, et projeter les points du semis (le dataframe arbres_intramuros) dans les cellules crƩƩes (cāest Ć©galement une faƧon dāapprocher une carte de densitĆ© 2D de faƧon discrĆØte)
Afin de ne pas faire une grille qui couvre la boĆ®te englobante de la zone dāĆ©tude, mais que la grille ne couvre que les gĆ©omĆ©tries de chaque quartier, on rĆ©alise une intersection entre la grille raster et lāenveloppe des quartiers.
N.B. ce nāest pas la faƧon canonique de rĆ©aliser un raster, pour cela il faut utiliser la library raster et rĆ©aliser le raster Ć partir lāaide de la fonction
rast1 <- st_make_grid(quartiers, square = TRUE, n=40, what="polygons") ## raster de 1600 cellules
plot(quartiers$geometry)
plot(rast1, add=T)
#restriction du raster à la geométrie des quartiers par intersection avec son enveloppe
envelop_qu <- st_union(quartiers)
plot(envelop_qu)
rast2 <- st_intersection(envelop_qu, rast1) ## raster de 1065 cellules
# on retransforme l'objet rast2 qui est une collection (sfc) en objet sf
rast2 <- st_sf(rast2)
plot(rast2)
#on récupère les points contenus dans les cellules
predicat_intersect <- st_contains(rast2,arbres_intramuros)
#on compte le nombre de points par cellules avec la fonction length
rast2$nb_arbres <- unlist(lapply(predicat_intersect, length))
hist(rast2$nb_arbres)
plot(rast2)
A quoi correspondent les cellules jaunes ?
On peut rƩitƩrer le processus en changeant la taille de la grille ou le type de grille (hexagonale)
rast1 <- st_make_grid(quartiers, square = TRUE, n=80)
rast2 <- st_intersection(envelop_qu, rast1)
rast2 <- st_sf(rast2)
predicat_intersect <- st_contains(rast2,arbres_intramuros)
rast2$nb_arbres <- unlist(lapply(predicat_intersect, length))
plot(quartiers$geometry, border="white", bgc="#222222" )
plot(rast2, border=NA, key.pos=4, add=T)
plot(quartiers$geometry, border="white", bgc="#222222", add=T, alpha=0.8,lwd=0.3,key.pos=4)
rasthex <- st_make_grid(quartiers, square = FALSE, n=60, what="polygons")
rasthex <- st_sf(rasthex)
#les hexagones s'arrètent à l'evloppe, pas besoin de la créer
predicat_intersect <- st_contains(rasthex,arbres_intramuros)
rasthex$nb_arbres <- unlist(lapply(predicat_intersect, length))
plot(quartiers$geometry,bgc="#222222")
plot(rasthex, key.pos=4, add=T, lwd=0.2)
plot(quartiers$geometry, border="white", bgc="#222222", add=T, alpha=0.8,lwd=0.3,key.pos=4)
On peut noter que le changement de dĆ©coupage de lāespace dĆ» Ć la taille ou la forme des cellules change lāaspect de la carte.
Pour calculer la densitĆ© dāarbres par quartier, nous devons compter le nombre dāarbres par quartier, et diviser ce nombre par la surface du qartier.
Nous avons dĆ©jĆ ajoutĆ© un attribut nb_arbres Ć lāobjet quartiers. Lāobjet quartiers contient aussi un attribut surface, issu des donnĆ©es brutes. On doit dāabord vĆ©rifier que cette valeur est correcte (la projection des donnĆ©es est Ć©quivalente) en calculant lāaires des quartiers et en les comparant Ć lāattribut prĆ©-existant
ecarts <- as.numeric(st_area(quartiers)) - quartiers$surface
ecarts
## [1] 1.877930e-05 8.108851e-06 6.117160e-06 5.210517e-06 4.029163e-06
## [6] 5.435315e-06 6.035494e-06 6.569724e-06 7.249997e-06 5.676586e-06
## [11] 8.165312e-06 4.341156e-06 6.586371e-06 9.512180e-06 1.055695e-05
## [16] 8.231029e-06 1.244282e-05 1.792773e-05 1.491560e-05 9.088253e-06
## [21] 6.806862e-06 1.615961e-05 1.833285e-05 5.744630e-06 1.758547e-05
## [26] 2.283882e-05 1.835986e-05 2.939277e-05 2.566329e-05 1.731119e-05
## [31] 1.683447e-05 2.600718e-05 1.549569e-05 1.163094e-05 8.670206e-06
## [36] 1.081359e-05 2.129877e-05 9.055191e-06 1.324248e-05 1.921703e-05
## [41] 1.531083e-05 1.852063e-05 2.501090e-05 2.072949e-05 1.291232e-04
## [46] 1.569642e-04 4.206994e-05 2.605096e-05 2.508820e-05 6.487360e-05
## [51] 4.781876e-05 1.511956e-05 2.415944e-05 2.820021e-05 2.937065e-05
## [56] 3.764336e-05 6.029010e-05 3.513182e-05 3.172900e-05 5.560508e-05
## [61] 1.362823e-04 1.181327e-04 6.534485e-05 3.181025e-05 3.170897e-05
## [66] 2.902956e-05 3.089476e-05 2.964167e-05 4.153443e-05 3.538677e-05
## [71] 2.405327e-05 2.895016e-05 2.844469e-05 5.151751e-05 3.934582e-05
## [76] 2.758135e-05 1.779071e-05 3.266637e-05 3.387243e-05 4.496775e-05
Les Ć©carts sont trĆØs faibles (de lāordre de \(10^{-5}m^2\)) , suffisamment pour que lāattribut surface soit directement utilisĆ© pour le calcul de la densitĆ©.
On crĆ©e une variable dens dans lāobjet quartier dont la valeur est le ratio entre la variable surface et la variable nb_arbres. On aura vĆ©rfiĆ© auparavant quāil nāy a pas de valeurs manquantes dans lāobjet , avec la fonction anyNA()
anyNA(quartiers)
## [1] FALSE
quartiers$dens <- quartiers$nb_arbres / quartiers$surface
plot(quartiers["dens"], main="DensitƩ d'arbres par quartier")
Connaissant lāimplantation des arbres et la formes des quartiers , que dire des densitĆ© des quartiers sud-ouest et sud-est ?
Que nous dit la comparaison entre la carte choroplète de nombre et celle de densité ?
Quel est le quartier le plus dense en arbres ?
leplusdense <- top_n(quartiers,1,dens)
leplusdense$l_qu
## [1] "PĆØre-Lachaise"
Pour utiliser des symboles proportionnels , nous devons utiliser un package spĆ©cifique : cartography. Dans ce package , lāajout de couche au dessin courant est activĆ© par dĆ©faut (add=TRUE)
library(cartography)
plot(quartiers$geometry, bgc="#888888", border="white", lwd=0.3)
propSymbolsLayer(quartiers,var = "nb_arbres", inches = 0.15, legend.pos = NA)
# couches de disposition des lƩgendes et du texte
layoutLayer(title = "Nombre d'arbres Ć Paris par quartier",
frame =F, north = TRUE, author="M2 IGAST2019-2020",
sources = "Opendata.paris.fr 2019",
scale = 50)
legendCirclesSymbols(pos = "topleft", title.txt = "nombre d'arbres", title.cex = 0.8, cex = 1,
border = "black", lwd = 1, values.cex = 0.6, var=c(50,1000,3500,7000), inches=0.15,
col = "#E84923", frame = FALSE, values.rnd = 0, style = "e")
Le package cartography est un excellent package, que je recommande vivement.